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. Quantum tomography is a procedure to determine the quantum state 

■~>^ ' of a physical system, or equivalently, to estimate the expectation value 

, of any operator. It consists in appropriately averaging the outcomes of 

the measurement results of different observables, obtained on identical 
copies of the same system. Alternatively, it consists in maximizing an 
appropriate likelihood function defined on the same data. The procedure 
can be also used to completely characterize an unknown apparatus. Here 
we focus on the electromagnetic field, where the tomographic observables 
are obtained from homodyne detection. 



■ Tomography. 
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1. Introduction 

The properties of each physical system are, by definition, completely deter- 
mined by its quantum state. Its mathematical description is given in form 
of a density operator g. Bohr's principle of complementaritj "'^ , which is in 
many ways connected with the uncertainty relationJ^, forbids one to recover 
the quantum state from a single physical system. In fact, the precise knowl- 
edge of one property of the system implies that the measurement outcomes 
of the complementary observables are all equiprobable: the properties of 
a single system related to complementary observables are simultaneously 
unknowable. Moreover, the no-cloning principle ' precludes to obtain many 
copies of a state starting from a single one, unless it is already known. Hence, 
complementarity and no-cloning prevent one to recover a complete infor- 
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mation starting from a single quantum system, i.e. to recover its state. The 
only possibility is to recover it from multiple copies of the system. [Notice 
that, if the multiple copies are not all in the same quantum state, we will re- 
cover the mixed state of the ensemble] . Given N copies of a system, we can 
either perform a collective measurement on all (or on subsets), or perform 
measurements separately on each system and combine the measurement 
results at the data analysis stage. Even though the former strategy would 
probably increase the speed of the statistical convergence of the measured 
state to the true one, it is quite impractical. Tomography thus adopts the 
latter strategy, which is the simplest to perform experimentally. 

What is quantum tomography? It is the name under which all state re- 
construction techniques are denoted. It derives from the fact that the first 
tomographic method (see Sec. Ujl employed the same concepts of Radon- 
transform inversion we find in conventional medical tomographic imaging. 
Since then, better methods have evolved which eliminate the bias that the 
Radon-transform necessarily entails. These fall into two main categories: the 
plain averaging method and the maximum likelihood method. As will be 
seen in detail, the first method requires a simple averaging of a function cal- 
culated on the N measurement outcomes Xn of the homodyne quadratures 
X^^. Thus, the statistical error which affects the estimated quantity can 
be easily evaluated through the variance of the data. The second method, 
i.e. the maximum likelihood method, is based on the assumption that the 
data we obtained is the most probable. Hence, we need to search for the 
state that maximizes the probability of such data, i.e. the state g for which 
YiiLi ct>n{^n\Q\xn)<pri is maximum, where <p„{xn\Q\xn)ct>„ is the probability 
of obtaining the result Xn when measuring the quadrature X^,^ (which has 
eigenstates la;)^^). 

Their involved mathematical derivation has given these tomographic 
techniques a false aura of being complicated procedures. This is totally 
unjustified: the reader only interested in applying the method can simply 
skip all the mathematical details and proceed to Sec. where we present 
only the end result, i.e. the procedure needed in practice for a tomography 
experiment (the experimental setup is, instead, given in Sec. I2.1|l . 

The chapter starts by introducing the method of homodyne tomography 
in Sec. [3 along with the description of homodyne detectors, noise deconvo- 
lution and adaptive techniques to reduce statistical errors. Then, in Sec. |31 
we present the Monte Carlo integration methods and the statistical error 
calculations that are necessary for the plain averaging technique. In Sec.^ 
the maximum likelihood methods are presented and analyzed. In Sec. El 
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the stcp-by-step procedure to perform in practice a tomography experi- 
ment is presented. In Sec. El a tomographic method to calibrate (i.e. com- 
pletely characterize) an unknown measurement device is presented. Finally, 
in Sec.[3 a historical excursus on the development of quantum tomography 
is briefly given. 

2. Homodyne tomography 

The method of homodyne tomography is a direct application of the fact that 
the displacements operators = ^aa'' -a a ^^^^ complete orthonormal 

set for the linear space of operators. Recalling that the scalar product in 
a space of operators takes the Hilbert-Schmidt form {A\B) —Tr[A^B], this 
means that 

f — TT[AV^ia)]V{a)^ /^"dr^Tri^ e^'-^^le"'''^* ,(1) 

JC I" Jo J-oo 4 

where the polar variables a = —ir e"^/2 were used in the second equal- 
ity. Upon introducing the probability p{x,(j)) = (p{x\g\x)(i, of obtaining x 
when measuring the quadrature = (a^'e''^ + ae"*"*)/2, one obtains the 
tomographic formula 

(A) = Tr[Ag] = dx p{x, <j>) Ka{x, </.) , (2) 

Jo J-oo 

where 

/+°° Irl 
dr^-^TT[A e-(^*~-)] , (3) 
-oo ^ 

defines the kernel of homodyne tomography. In the case of the density 
matrix reconstruction in the Fock basis \n) (i.e. when A — \n){m\), the 
kernel function i^ 

i^^(x,0) = 2e'(™-")^He-^ y tllf ^ ) (4) 

x{2j+n-m + iy. Re[(-l)"-"I?_2(2,+„-™+2)(-2^a;)] , 

where Re denotes the real part and T>i{x) denotes the parabolic cylinder 
function (which can be easily calculated through its recursion formulas). 

The multimode case is immediately obtained by observing that the 
quadrature operators for different modes commute, so that for an oper- 
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ator Am (acting on the Hilbert space of M modes) we find 




xii'AM(a;i, 01, ■ • ■ , xm, 4>m) , 



(5) 



wfiere p{xi, - ■ ■ ,xm,4'm) the joint probability of obtaining the results 
{xjn} when measuring the quadratures {X^^}, and where 



KA,A^u<t'i,---)= / dn---drM n ^Tr[AM e^'""^^^"^"'")] . (6) 



However, such a simple generalization to multimode fields requires a sepa- 
rate homodyne detector for each mode, which is unfeasible when the modes 
of the field are not spatio-temporally separated. This is the case, for exam- 
ple of pulsed fields, for which a general multimode tomographic method 
is especially needed, because of the problem of mode matching between 
the local oscillator and the detected fields (determined by their relative 
spatio-temporal overlap), which produces a dramatic reduction of the over- 
all quantum efficiency. A general method for multimode homodyne tomog- 
raphy can be found^I that uses a single local oscillator that randomly scans 
all possible linear combinations of incident modes. 

2.1. Homodyne Detection 

The balanced homodyne detectoJ^ measures the quadratures = (a^e^'^-f 
ae~^'^)/2. The experimental setup is described in Fig.^ The input-output 
transformations of the modes a and b that impinge into a 50-50 beam- 
splitter are c = (a + b)/ \/2, d = (a — b)/ \/2 where c and d are the two beam- 
splitter output modes, each of which impinge into a different photodetector. 
The difference of the two photocurrents is the homodyne detector's output, 
and thus is proportional to c^c — d'^ d = a^b + b^'a. In the strong local 
oscillator limit, with mode b in an excited coherent state (|/3| ^ 1), the 
expectation value of the output is Ih oc {a^)f]+ {a)(i* which is proportional 
to the expectation value of the quadrature X^, with <f) the relative phase of 
the local oscillator. 

A detector with non-unit quantum efficiency is equivalent^ to a per- 
fect rj — 100% detector, preceded by a beam-splitter with transmissivity 
rj. Inserting two beam-splitters in front of the two photodiodes of the ho- 
modyne scheme, the modes c and d evolve as c' = y/r] c + y/1 — rj u and 
d' = yff\ d -\- \J\ — T] V, where u and v are vacuum noise modes. The homo- 
dyne output, is now proportional to c'''c' — d'^'d' , i.e. to L = r] (a^6 + b^a) + 
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Fig. 1. Homodyne detector. The input signal (in mode a) is mixed by a 50-50 beam- 
splitter (BS) with a strong local oscillator (LO), which is coherent with the input field 
and is in a strong coherent state. The relative phase (f> between the signal and the LO 
must be known and should be varied in [0, tt] with uniform probability. Two identical 
high efficiency linear photodetectors Pi and P2 measure the field. The photocurrents are 
then accurately subtracted electronically yielding the output Ih ■ Since the LO amplifies 
the weak quantum signals of the input, one can use high efficiency detectors that work 
only with strong signals. 

{l-ri){u''u-v^v) + ^{l - f])f] /2[a{u'< -v'^)+b{u^ +v'')+a'^ {u~v)+b^u+v)]. 
As before, we take the limit |/3| ^ 1 of strong pump in b, and rescale the 
output difference photocurrent by 2|/3|?7, obtaining 



where the modes u and v are in the vacuum state. Since the quadrature 
outcome for each vacuum state is Gaussian-distributed with variance 1/4, 
this means that the distribution of the noisy data are a convolution of the 
clean data with a Gaussian of variance = (1 — i])/{4ri), namely 



2.2. Noise deconvolution 

The data-analysis procedure can be modified to yield the result we would 
obtain from perfect detectors, even though the data was collected with noisy 
oneJ^. In fact, depending on which operator A we consider and on the value 
of the quantum efficiency rj, the noise may be numerically deconvolved. The 
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output of the noisy homodyne is distributed according to Eq. (jS)) , and one 
can rewrite Eq. |2l as follows 

{A) = r ^ r dx p^{x,<t>) r dr M e^-'^'/^ Tr[v4 e'-(^*--)] , (9) 

>/ — CXD o — OO 

where Priix, 4>) is the probability of the noisy data. In the case when all the 
integrals are convergent, the noise inversion can be performed successfully. 

It is clear the possibility of noise deconvolution depends on the quantum 
efficiency of the detectors and the operator to be estimated. For example, 
there is a bound r] > 50% for the reconstruction of the density matrix in 
the Pock basis (i.e. for A = \n){m\). In fact, one can see that for rj < 50% 
Eq. @ has an unbounded kernel. Notice that actual homodyne detectors 
have efficiencies ranging between 70% and 90%. 

2.3. Adaptive tomography 

Adaptive tomographjl^ exploits the existence of null estimators to reduce 
statistical errors. In fact, the addition of a null estimator in the ideal case 
of infinite statistics does not change the average of the data since, by def- 
inition, the mean value of a null estimator is zero. However, it can change 
the variance of the data. Thus, one can look for a procedure to reduce the 
variance by adding suitable null functions. 

In homodyne tomography null estimators are obtained as linear combi- 
nations of the following operators 

AfkAX^) - ^±^{k+2+2n)v ^ k,n>0. (10) 

One can easily check that such functions have zero average over ip, indepen- 
dently on g. Hence, for every operator A one actually has an equivalence 
class of infinitely many unbiased estimators, which differ by a linear combi- 
nation of functions Afk,n{Xip). It is then possible to minimize the rms error 
in the equivalence class by the least-squares method. This yields an opti- 
mal estimator that is adapted to the particular set of experimental data. 
Examples of simulations of the adaptive technique that efficiently reduce 
statistical noise of homodyne tomographic reconstructions can be found in 
Ref.El. 

3. Monte Carlo methods for tomography 

In this section we will very briefly review the basics of the Monte Carlo 
integration techniques that are needed and we show how to evaluate the 
statistical error bars of the tomographically estimated quantities. 
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A tomographic technique is based on an integral of the form 

r+oo 

F= / dxpix)f{x) , (11) 



where p{x) is a probability. Since we have experimental outcomes {xn, n — 
1, • • • N} distributed according to the probability p{x), we sample the inte- 
gral ((nT) using 



+00 2 



N 

dxp{x)f{x)^ lun (12) 

n=l 

For finite N, the sum will be an unbiased estimator for the integral, af- 
fected by statistical errors only (which can be made arbitrarily small by 
increasing N). The central limit theorem guarantees that the finite sum 
Fjv = J2n=i f{^n)/N is a statistical variable distributed as a Gaussian (for 
sufficiently high N) with mean value F and variance 



1 ^ 1 2/ T?\ 



n=l j=l 

Hence, the tomographic estimated quantity converges with a statistical er- 
ror that decreases as 1/ VTV- It can be estimated from the data as 

1 ^ 

s'{Fn) = j^^Y.^F^-^f- (14) 

[Remember that the factor iV — 1 in the variance denominator arises from 
the fact that we are using the experimental estimated mean value m in 
place of the real one F .] The variance of the statistical variable 'mean m' 
is then given by a'^(m) = a'^{Fjq)/N, and thus the error bar on the mean 
m estimated from the data is given by 



N ' N{N- 

From the Gaussian integral one recovers the usual statistical interpretation 
to the obtained results: the "real" value F is to be found in the interval 
[m — e,m + e] with ~ 68% probability, in the interval [m — 2e,m + 2e] with 
~ 95% probability and in [m ~ 3e,m + 3e] with unit probability. 

In order to test that the confidence intervals are estimated correctly and 
that errors in the data analysis or systematic errors in the experimental data 
do not undermine the final result, one may check the F^ distribution, to 
see if it actually is a Gaussian distribution. This can be done by comparing 
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a histogram of the data to a Gaussian, or by using the test. Notice that 
when we have very low statistics it may be useful to use also bootstrapping 
techniques to calculate the variance of the data. 

For a more rigorous treatment of the statistical properties of q uantum 
tomography, and also some open statistical questions, see Ref. 

4. MzLximum likelihood tomography 

The maximum likelihood tomography is based on the as sum ption that the 
data obtained from the measurements is the most likeljUJJ. In contrast to 
the plain averaging method presented above, the outcome is not a simple 
average of functions of the data, but a Lagrange-multiplier maximization is 
usually involved. The additional complexity introduced is compensated by 
the fact that the results are statistically less noisy. Estimation of operator 
expectation values is, however, indirect: one must first estimate the state g 
and then calculate the expectation value as Tr[pyl]. 

Consider a known probability distribution P'y{x) parametrized by a pa- 
rameter 7 (which may also be a multidimensional parameter). We want to 
estimate the value of 7 from the data set {xi, • • • , xn}. The joint probabil- 
ity of obtaining such data is given by the likelihood function 



The maximum likelihood procedure consists essentially in finding the 
7o(a:i, • • • , Xn) which maximizes the likelihood function C{xi, • • • , xn', 7). 
Equivalently, it may be convenient to maximize its logarithm 
log£(xi, • • • , xat; 7), in order to convert into a sum the product in Eq. H16|l . 
Usually, various constraints are known on the parameters 7, which can be 
taken into account by performing a constrained maximization. The confi- 
dence interval for the estimated 70 can be evaluated from the data using a 
bootstrapping technique: we can extract a rough estimate of the probability 
distribution of the {xi} from the data set, generate M simulated sets of N 
data points, and repeat the procedure to obtain a set of M parameters 7o™^ • 
Their variance estimates the variance of the reconstruction. Moreover, if a 
sufficiently large data set is present, we can attain the Cramer- Rao bound 
cr^ ^ 1/iVF-y, where is the Fisher information relative to Pj{x), i.e. 



£{xi 



(16) 




(17) 
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Since the Cramer-Rao bound is achieved only for the optimal estimatoi^l 
the maximum likelihood is among the best (i.e. least statistically noisy) 
estimation procedures. 

Th e rn aximum likelihood method can be extended to the quantum 
domair^^. The probability distribution of a measurement is given by the 
Born rule as pi =Ti[Ilig] where {Hi} is the positive operator- valued mea- 
surement (POVM) that describes the measurement. Thus we need to max- 
imize the log-likelihood function L(g) = log Tr[ni£i] over the set of den- 
sity operators g. In the case of finite Hilbert space, L{g) is a concave function 
defined on a convex set of density operators: its maximum is achieved on 
a single point or on a convex subset. The main difficulty of this procedure 
consists in finding a simple parameterization for the density matrix, that 
enforces both the positivity and the normalization Tr[p] = 1. The former 
is guaranteed by requiring that g ~ T'^T, the latter must be taken into 
account through an appropriate Lagrange multiplier. In order to employ 
the minimum number of parameters, it is sufficient to consider T as an 
upper complex triangular matrix with nonnegative diagonal elements — so 
called Cholesky decomposition. This decomposition achieves minimal pa- 
rameterization (up to the normalization condition), as it requires (P real 
parameters for a d x d Hermitian matrix. Thus, in practice we need to 
maximize the operator Lx[g] = logTrpiT^T] — A[r^r], where A is a 
Lagrange multiplier that accounts for the normalization. By expressing g 
in terms of its eigenstates as ^3 = X)m ymlV'm)(V'm|, the condition for the 
maximum, dL\/dym = 0, becomes 

J2{ym{^mm^m)/Tr[gn,]}~ Xy„,^0 Vm. (18) 

i 

Multiplying both members by and summing over m, through the Born 
rule and the normalization of g, we find that A is equal to the number 
of measurements employed. Thus, we are left with the problem of finding 
the maximum of the d^-parameter function La=jv[£' — T^T], which can 
be tackled with conventional nu mer ical techniques such as expectation- 
maximization or downhill simple?!^. By using the ML method only small 
samples of data are required for a precise determination, even in the pres- 
ence of low quantum efficiency at the detectors. However, we want to em- 
phasize that such method is not always the optimal solution of the tomo- 
graphic problem, since it suffers from some major limitations. Besides being 
biased due to the Hilbert space truncation — even though the bias can be 
very small if, from other methods, we know where to truncate — it cannot 
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be generalized to the estimation of any ensemble average, but just of a set 
of parameters from which the density matrix depends. In addition, for the 
multi-mode case, the method has exponential complexity versus the number 
of modes. 

5. Tomography for dummies 

In this section we just give the step-by-step procedure to implement a to- 
mography experiment, employing all the results obtained in the previous 
sections. 

• Plain averaging method: 

(1) Calculate the Kernel function Ka for the operator A whose expectation 
value we want to estimate through Eq. 0. For example, to estimate 
the density matrix in the Fock basis, we need the Ka defined in Eq. Q). 

(2) The experimental apparatus, described in Sec. 12. II yields a set of N data 
points : each datum is composed by the quadrature phase (/)„ 
that was measured and by the corresponding measurement result Xn ■ 

(3) Evaluate ^„ X^(x„, (/)„). In the limit N ^ oo this average yields 
the expectation value (A) we are looking for. 

(4) For finite N, we can estimate the purely statistical error on the result 
through Eq. H15|l , replacing m with the average obtained at the previous 
point and Fn with the nth Kernel function evaluation, KA[xn,X^^)]. 

Further data massaging is also possible: we can employ adaptive tomogra- 
phy to reduce the statistical noise (see Sec I2.3|l . Moreover, we can remove 
the detector noise due to homodyne measurements with non unit quantum 
efficiency ry, as long as i] > 1/2 (see Sec. l2.2() . 

• Maximum likelihood method: 

(1) Parametrize the unknown quantum state through the upper triangular 
dxd matrix T as g = T^^T. 

(2) Use the same experimental apparatus (homodyne detection) to obtain 
N data points {(/)„, Calculate the log likelihood function on the 
experimental data as logX!^=i (l>J\^n\T^T\xn) 4>ri- 

(3) Numerically maximize this quantity over the d? parameters of T with 
the additional constrain Ty[T^T] = 1. This maximum is achieved on 
our best estimate for the state g — T^T. 

(4) The confidence intervals for our estimation can be obtained using boot- 
strapping techniques, or employing the Cramer-Rao bound of Eq. IjlTfl. 
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6. Quantum calibration of measurement devices 

In this section we review the method to measure the POVM of an unknown 
measurement apparatus presented in Ref. The method is based on an- 
alyzing the correlations in measurements on a bipartite system: one of the 
two parts is fed into the unknown apparatus A, while the other is measured 
with a known set B of detectors that measures a quorum of observables 
(see Fig. 121). As will be shown in the following, there is ample freedom in 
the choice of both the input bipartite states and the set of observables. The 
procedure is repeated many times and the joint measurement outcomes are 
analyzed using the tomographic algorithms described above, which (in the 
limit of infinite input data) yield the POVM of the unknown apparatus. 
For finite data, the reconstructed POVM will be affected only by statis- 
tical errors which can be easily estimated. For the sake of illustration, a 
Monte-Carlo simulation of the procedure is given at the end of this section. 
It aptly illustrates the advantage of using maximum likelihood techniques 
over plain averaging: the maximum likelihood reconstruction is significantly 
less noisy. 




Fig. 2. (Left) Experimental setup to determine the POVM of the unknown measure- 
ment apparatus A: one part of the bipartite input state R is sent to the apparatus A 
which yields the measurement result n; the other part (with quantum state g„) is sent 
to the known detector B which performs a projective measurement of an observable B^. 
from the complete set {-Bfe} yielding the result m{k). The joint measurement results are 
processed using a tomographic algorithm to obtain the POVM {n„} of A. (Right) Exam- 
ple of application of the scheme to the radiation field. The bipartite state R is generated 
via a non-linear crystal through spontaneous parametric down-conversion. The tomog- 
rapher B is, in this case, a homodyne detector (HD) which measures the quadratures, a 
complete set of observables. 

The following simple example illustrates how the procedure works. Sup- 
pose we want to evaluate the POVM of a von Neumann measurement of 
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the observable O which acts on a d-dimensional Hilbert space Ha and has 
spectral decomposition X^n maximally entan- 
gled input state l^*) = J2i=i N)N)/v^: which lives in the space Ha <8) Hb- 
In fact, this state can be also written as 




where * denotes complex conjugation with respect to the basis \i). It is 
obvious from Eq. H19|) that the outcome o„ at detector A (corresponding 
to the state |o„) in Ha) means that the state gi„ = |o*)(o*| in Hb im- 
pinges in detector B. The POVM can be recovered using tomographical 
state reconstruction at B, since in this simple case n„ = gi* . 

It is not difficult to generalize the above example to arbitrary POVMs 
and measurement procedures. Let the unknown apparatus A be described 
by the POVM {n„} we want to estimate, and let the apparatus B mea- 
sure the quorum observables Ok described by the von Neumann projections 
{\km) {km\} (with basis for all k). From the Born statistical formula 

we can derive the state that impinges into the known detector B if the un- 
known detector A gave result n for the measurement on the initial bipartite 
state R, as 

_ Tri[(n„^IL)i?] 
^" - Tr[(n„ ® t)R] ■ ^^"^ 

It describes the state reduction at B stemming from a measurement at A 
with outcome n. The denominator is the probability p{n) of obtaining the 
result n at B. The state gi„ contains some information on the POVM element 
n„. It can be recovered by introducing the map Tl{X) = Tri[(X(g) so 
that Eq. rewrites as Qn = TZ[Un/p{n)]. This implies that the POVM 
can be recovered as n„ = p{n)TZ~^ {gn) , where the map TZ depends only on 
the input state R: the input state R allows the POVM reconstruction if the 
inverse map TZ~^ exists. This condition can be cast in a more transparent 
form by rewriting the map 7?, in a m ultip licative form via isomorphism 
between operators oti H®H and map»^. We can obtain an operator of 
this form by considering S ~ R^^ , i.e. the partial transposition on the first 
space of the input state R. In fact, taking two operators X and Y such that 
Y = Ti{X), we see that 

Yu = ^X,,(^|7^(|J)(fc|)|0 =J2Xjk{R^')jk.u , (21) 

jk jk 
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where Yu = {i\Y\l), X,k = {j\X\k}, and {R^')jk,ii = {j\{i\R'^'\k)\l), the 



set {|n)} being a basis in Ti.. In matrix notation (considering jk and il as 
collective indexes), Eq. H21|l rewrites as Y = SX. It follows immediately 
that the map TZ is invcrtiblc if S^^ exists so that X — S^^Y. In this case 
we say that the input state R is faithfuM Since invertibility is a condition 
satisfied by a dense set of operators, the set of input states R that allow the 
POVM reconstruction is also dense, i.e. almost any bipartite state will do. 
In particular, all Gaussian bipartite states — with the trivial exception of 
product states — are faithful To recapitulate: in order to check whether 
the state £»„ allows to obtain the POVM (i.e. whether the input state R is 
faithful) we must verify that the operator (i?"^i )jk,ii is invertible when jk 
and il are considered as collective indexes. As an illustration of this check, 
take the simple example given above: the state |^) — jyfd is faithful 

since |^)(\E'|"^i — ^Yliij is invertible: it is a multiple of the swap 

operator E = Y,^ 

To recover from the measurements at B (and hence the POVM if 
the input R is faithful), we can use the quantum tomographic techniques 
described in the previous sections. If we employ the plain averaging tech- 
nique, we may recover the density matrix elements Qij in some basis and 
then calculate the POVM using the inverse map 7?."^, as 



where the inverse of must be calculated considering jfc and il as col- 
lective indexes. On the other hand, if we employ maximum likelihood we 
may directly maximize the probability of acquiring the data we obtained 
from the measurements-'^ i.e. the joint probability pk{n,m) =Tr[(n„ (g) 
\km){km\)R\- Equivalently, one can maximize the logarithm of this quan- 
tity and consider simultaneously all the N joint measurement outcomes 
{ni, TOi}, • • • , {nNiiTiN} of the quorum operators O^ci) at detector A and 
of the unknown detector B. Thus, the POVM {n„} is the one that maxi- 
mizes the quantity 



with the additional constraints n„ ^ and ^„ n„ = 1. Other prior knowl- 
edge on the quantities to be estimated can be easily introduced adding 
further constraints to the maximization. Also in this case it is possible to 
take into account a known source of noise at the detector B: if we replace 



{3\n^\k)^p{n)J2g^'> (i?^0-,!.^ 



il 



N 




(22) 
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the term \kmi){kmi \ in Eq. with the noise-evolved A/'(|fcmi)(fcm,l)7 then 
the maximization yields the POVM that maximizes the noisy measurement 
results. 

For the sake of illustration, we give a Monte-Carlo simulation of the 
calibration procedure in which we recover the POVM of a simple inefficient 
photodetectoi-'^'^. An inefficient photodetector is aptly modeled by a perfect 
photodetector (which is a device which measures the observable "number 
of photons" a^a = ""I"") (""D' preceded by a beam-splitter with a trans- 
missivity equal to the quantum efficiency rj of the detector. Possible dark 
counts can be considered by feeding the other beam-splitter port with a 
thermal state with n average photons. In this case, the theoretical POVM 
is given by 

oo 

n„ = ^|p)(p| (23) 

CXD min(p,fc+n) / s / i \ / 7 , \ 

fe=o i=o V-^/V ' / \J / 

Since this POVM is diagonal in the Fock basis, we can limit the reconstruc- 
tion to the diagonal elements. As input state R we employ a twin beam state 
\TB), i.e. the result of spontaneous parametric down-conversion: 

\TB) ^ v/T~^^r|m).|m), , (24) 

where ^ is the parametric amplifier gain and \m)a and \m)h are Fock states of 
the modes a and b that impinge in the detectors A and B respectively. This 
is a faithful state since \TB){TB\^^ = (1- |^|2)i;^°*'^(g)^*''*^ (where is the 
swap operator) is invertible. The photon counter measures the mode a at 
position A, while homodyne detection with quantum efhciency rjh measures 
the mode b at position B acting as tomographer (see Fig.|2J). Since only the 
diagonal part of the POVM is needed, we can use a homodyne detector with 
uniformly distributed local oscillator phase. [A phase-controlled homodyne 
detector would allow to recover also the off-diagonal elements of the POVM, 
ensuring a complete characterization of the device.] 

In Figs. 01 and ^ we present the results of the POVM reconstruction 
deriving from the two tomographic methods described above (simple av- 
eraging and maximum likelihood, respectively). The convergence of the 
maximum likelihood procedure is assured since the likelihood functional 
C is convex over the space of diagonal POVMs. However, the convergence 
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speed can become very slow: in the simulation of Fig.0|a mixture of sequen- 
tial quadratic programming (to perform the constrained maximization) and 
expectation-maximization techniques were employed. From the graphs it is 
evident that the maximum likelihood estimation is statistically more effi- 
cient since it needs much less experimental data than tomography. This is 
a general characteristic of this method, since if the optimal estimator (i.e. 
the one achieving the Cramer-Rao bound) exists, then it is equal to the 
maximum likelihood estimator An added bonus, evident from Eq. H22|l . 
is that the maximum likelihood recovers all the POVM elements at the 
same time additionally increasing the statistical efficiency. On the other 
hand, the tomographic reconstruction is completely unbiased: no previous 
information on the quantity to be recovered is introduced. 

This simulated experiment uses reali stic parameters and is feasible in 
the lab with currently available technologjU^. The major experimental chal- 
lenge lies in the phase matching of the detectors, i.e. in ensuring that the 
modes detected at A and B actually correspond to the modes a and h of 
the state \TB). 



7. History of quantum tomography 

In this section a brief historical perspective (see alsc^^l) on quantum 
tomography is presented. Already in 1957 FancPI 

stated the problem of 

quantum state measurement, followed by rather extensive theoretical work. 
It was only with the proposal by Vogel and Risken however, that homo- 
dyne tomography was born. The first experiments followed^ by showing 
reconstructions of coherent and squeezed states. The main idea at the basis 
of these works, is that it is possible to extend to the quantum domain the 
algorithms that are conventionally used in medical tomographic imaging 
to recover two-dimensional distributions (say of mass) from unidimensional 
projections in different directions. However, these first tomographic meth- 
ods are unreliable for the measurement of unknown quantum states, since 
some arbitrary smoothing parameters have to be introduced. 

A new approach to optical tomography was then proposed^^SSI which 
allows to recover the quantum state of the field g (and also the mean values 
of system operators) directly from the data, abolishing all the sources of 
systematic errors. Only statistical errors (that can be reduced arbitrarily by 
collecting more experimental data) are left. Quantum tomography has b een 
then generalized to the estimation of arbitrary observable of the field^, 
to any number of modes'', and to arbitrary quantum systems via group 
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Fig. 3. (Above left) Theoretical value of the diagonals of the POVM elements (m|nn |m) 
of the inefficient photodetector described by Eq. 1231 . with parameters n = 1, rj = 80%. 
(Above right) Simulated reconstruction of the same quantity. The data are simulated as 
coming from an input twin-beam state \TB) with § = 0.88, and as being detected from 
a phase insensitive homodyne detector with quantum efficiency r/f^ = 90%. Here 5 X 10^ 
simulated homodyne measurements are employed. (Below) The same data is plotted 
separately for each POVM element to emphasize the error bars. They are obtained from 
the root-mean-square of the recovered POVM matrix elements. (The theoretical value is 
plotted as the thick dashed line.) Plain tomographic averaging with noise deconvolution 
has been employed here, since the noise map of inefficient homodyne detection can be 
inverted for r]f^ > 50%. 



theorjES with further improvements such as noise deconvolutiorPI adap- 
tive tomographic methods^, and the use of max-Hkchhood strategieJIH 
which has made possible to reduce dramaticaUy the number of experimen- 
tal data, with negligible bias for most practical cases of interest. The latest 
developments are based on a general method^, where the tomographic 
reconstruction is based on the existence of spanning sets of operators, of 
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Fig. 4. Maximum likelihood reconstruction of the same POVM of Fig. 13 with the same 
parameters, but here only 5 X 10^ simulated homodyne measurements are employed. 
The statistical error bars are obtained by bootstrapping, i.e. by calculating the variance 
using the data of 50 numerical experiments. Notice that the result is statistically less 
noisy than the results presented in Fig. U] even if here less measurements are employed: 
maximum likelihood is usually a better estimator. 

which group tomographjEHI ig j^st a special case. 
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